Load Library Packages

library(tidyverse)
library(sf)
library(mapview)
library(leafpop)
library(tigris)
library(DT)

Load Data

starbucksNC <- read_csv("https://raw.githubusercontent.com/libjohn/mapping-with-R/master/data/All_Starbucks_Locations_in_the_US_-_Map.csv") %>% 
  filter(State == "NC")
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   `Store Number` = col_double(),
##   `Facility ID` = col_double(),
##   `Food Region` = col_double(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
options(tigris_use_cache = TRUE)

durham_blocks <- tracts(state = "NC", county = "Durham", class = "sf", year = 2012)

durham_blocks <-  durham_blocks %>%  
  select(-6, -7, -8, -11, -12)

glimpse(durham_blocks)
## Observations: 60
## Variables: 8
## $ STATEFP  <chr> "37", "37", "37", "37", "37", "37", "37", "37", "37",...
## $ COUNTYFP <chr> "063", "063", "063", "063", "063", "063", "063", "063...
## $ TRACTCE  <chr> "001001", "000200", "000900", "000302", "000102", "00...
## $ GEOID    <chr> "37063001001", "37063000200", "37063000900", "3706300...
## $ NAME     <chr> "10.01", "2", "9", "3.02", "1.02", "10.02", "4.01", "...
## $ ALAND    <dbl> 2283723, 3021358, 1025539, 1722150, 4011880, 3348806,...
## $ AWATER   <dbl> 3268, 0, 0, 0, 3346, 0, 43208, 0, 121500, 0, 3875, 18...
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-78.88456 3..., MULTIPOL...

Wrangle

sbux_drm <- starbucksNC %>% 
  filter(City == "Durham") %>% 
  select(`Store Number`, `Ownership Type`,
         `Phone Number`, City, Zip, Longitude,
         Latitude)


# assign my coordinates to a degree system
# documenting my XY values in degrees
# use, for example 4326 (best for US) ; 4269 (for US Agency)
sbux_drm_sf <- st_as_sf(sbux_drm,
                        coords = c("Longitude", "Latitude"),
                        crs = 4326)
st_crs(starbucksNC)
## Coordinate Reference System: NA
st_crs(sbux_drm_sf)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(durham_blocks)
## Coordinate Reference System:
##   EPSG: 4269 
##   proj4string: "+proj=longlat +datum=NAD83 +no_defs"
# https://guides.library.duke.edu/r-geospatial/CRS
#       >    Projections / CRS > ESPG codes for CR's

Tranform to a Projected CRS

Use one of the NC Projections

# Transform to a Projection
durham_blocks <- st_transform(durham_blocks,
                              st_crs(2264))

sbux_drm_sf <- st_transform(sbux_drm_sf,
                            st_crs(durham_blocks)) 



st_crs(sbux_drm_sf)
## Coordinate Reference System:
##   EPSG: 2264 
##   proj4string: "+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"
st_crs(durham_blocks)
## Coordinate Reference System:
##   EPSG: 2264 
##   proj4string: "+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"

First Lookie

mapview(durham_blocks)
mapview(sbux_drm_sf)
mapview(list(durham_blocks, sbux_drm_sf))

join

shp_with_bux <- st_join(durham_blocks, sbux_drm_sf, join = st_contains, left = FALSE)

bux_within_shp <- st_join(sbux_drm_sf, durham_blocks, left = TRUE, join = st_within)
plot(st_geometry(shp_with_bux))

plot(st_geometry(bux_within_shp))

mapview(durham_blocks, col.regions = "grey",
        legend = FALSE,
        popup = popupTable(st_drop_geometry(durham_blocks), 
                           zcol = c(4,5),
                           row.numbers = FALSE,
                           feature.id = FALSE)) +
  mapview(shp_with_bux, 
          zcol = "NAME",
          layer.name = "Tract Name") +
  mapview(bux_within_shp,
          legend = FALSE,
          cex = 0,
          #cex = "AWATER",
          lwd = 3,
          alpha = .5,
          color = "red",
          zcol = "Store Number",
          popup =
            popupTable(st_drop_geometry(bux_within_shp),
                             zcol = c(1,2,3,5,10,11),
                           row.numbers = FALSE,
                           feature.id = FALSE))

Display Table

datatable(st_drop_geometry(bux_within_shp))
bux_dt <- st_drop_geometry(bux_within_shp) %>% 
  select(1, GEOID, TractName = NAME, Ownership = `Ownership Type`) %>% 
  group_by(TractName) %>% 
  summarise(TotalStores = n()) %>% 
  arrange(-TotalStores) 

datatable(bux_dt)

Choropleth

Shade Tract shape by number of Stores in the Tract

bux_tot_in_shape <- left_join(shp_with_bux, bux_dt, by = c("NAME" = "TractName"))

mapview(durham_blocks, col.regions = "grey",
        legend = FALSE,
        popup = popupTable(st_drop_geometry(durham_blocks), 
                           zcol = c(4,5),
                           row.numbers = FALSE,
                           feature.id = FALSE)) +
  mapview(bux_tot_in_shape, 
          zcol = "TotalStores",
          layer.name = "Total Stores by Tract",
          col.regions = c("royalblue1", "royalblue4")) +
  mapview(bux_within_shp,
          legend = FALSE,
          cex = 0,
          #cex = "AWATER",
          lwd = 4,
          color = "yellow",
          zcol = "Store Number",
          popup =
            popupTable(st_drop_geometry(bux_within_shp),
                             zcol = c(1,2,3,5,10,11),
                           row.numbers = FALSE,
                           feature.id = FALSE))